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Sub-Grid  Models  for  Large  Eddy  Simulation  of  Turbulent 
Combustion 


1.1  Introduction 

In  the  past  decade,  Large  Eddy  Simulation  (LES)  has  been  increasingly  and  successfully  applied 
to  both  premixed  and  non-premixed  reacting  flows  [1,  2,  3].  In  application  and  testing,  meth¬ 
ods  such  as  steady  and  unsteady  flamelet  modeling,  the  probability  density  function  (PDF) 
assumption,  and  level  set  tracking  have  been  shown  to  describe  combustion  in  a  LES  context 
correctly.  As  capable  as  such  methods  are,  however,  researchers  only  very  recently  have  begun 
to  apply  them  to  the  detailed,  multi-scale,  and  multi-physics  flows  that  arise  in  the  burners  and 
engines  found  in  modern,  industrially  relevant  equipment.  In  the  course  of  this  transition  of 
LES  from  a  scientifically  interesting  method  in  its  own  right  to  a  research  platform  useful  for 
studying  real  world  flows,  the  limits  of  current  methods  to  be  highlighted. 

More  specifically,  as  the  complexity  of  the  flows  being  considered  increases,  so  do  the  number 
of  involved  lengthscales,  the  details  of  the  relevant  chemistry,  and  the  extent  of  the  turbulence- 
chemistry  interactions.  In  LES,  where  the  large  scales  of  the  flow  are  resolved  while  the  small 
scales  are  modeled,  this  increased  complexity  means  that  more  accuracy  will  be  demanded  of 
sub-grid  models.  Sub-grid  models  for  momentum  have  been  studied  rigorously,  and  computa¬ 
tionally  efficient  methods  exist.  But  sub-grid  models  for  combustion-related  phenomena  have 
not  yet  demonstrated  the  ability  to  capture  the  physics  needed  to  describe  complex  flows  con¬ 
sistently  and  accurately.  For  example,  quantities  of  interest  such  as  flame  lift-off  heights  and 
heat  release  near  injection  nozzles,  often  are  predicted  incorrectly  using  state-of-the-art  codes. 
Furthermore,  NOx  and  soot  formation  processes,  which  occur  on  timescales  significantly  longer 
than  those  of  CO2  and  H2O  production,  cannot  be  captured  with  the  kind  of  steady  flamelet 
models  now  being  used  [4]. 

Part  I  of  this  project  was  proposed  in  response  to  these  challenges.  The  goal  of  part  I  has 
been,  broadly,  to  improve  the  predictive  capability  of  LES  of  turbulent  reacting  flows.  To  meet 
this  challenge,  more  specific  goals  were  set  for  each  of  the  three  combustion  regimes  of  industrial 
interest:  non-premixed,  premixed,  and  partially-premixed  flows. 

In  non-premixed  flows,  fuel  and  oxidizer  enter  a  combustion  chamber  separately.  These  com¬ 
ponents  are  combined  by  mixing  processes  to  sustain  chemical  reactions  within  the  chamber. 
In  turbulent  flows  where  small  eddies  can  influence  mixing  processes  strongly,  sub-grid  scale 
models  for  the  mixture  fraction  must  be  accurate  and  robust.  The  goal  of  the  non-premixed 
portion  of  this  project  was  to  assess  and  improve  models  for  the  scalar  variance  and  to  test  the 
limitations  of  these  models  when  used  in  complex  flame  simulations. 

In  premixed  flows,  reactions  most  often  occur  in  thin  flame  fronts.  These  reactions  induce 
front  propagation,  a  phenomenon  described  by  the  burning  velocity.  The  burning  velocity  is  a 
scalar  quantity  that  must  describe  the  effects  of  chemistry  as  a  function  of  flow  conditions.  In 
turbulent  flames,  where  eddies  can  influence  the  propagation  of  the  front,  the  burning  velocity 
also  must  describe  the  effects  of  turbulence  on  front  propagation.  Thus,  accurately  modeling 
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this  quantity  is  of  crucial  importance.  The  second  goal  of  part  I,  then,  was  the  development 
of  a  new  method  of  describing  front  propagation  in  premixed  turbulent  flows.  This  method 
should  be  able  to  account  for  variations  of  local  equivalence  ratio,  the  amount  of  unresolved 
local  turbulence,  and  the  local  resolution  of  the  numerical  flame  front  solution. 

Finally,  in  partially-premixed  flows,  chemical  reactions  can  occur  in  both  of  the  modes  de¬ 
scribed  above.  Regions  where  fuel  and  oxidizer  have  mixed  but  not  ignited  may  exist  locally 
within  the  flow,  such  as  in  strong  recirculation  zones  or  in  the  pre-burn  region  of  lifted  flames. 
Conversely,  purely  diffusive  burning  may  occur  in  regions  of  strong  shear  where  fuel  and  oxidizer 
streams  interact.  Such  flames  recently  have  received  the  attention  of  many  researchers  [5,  6], 
but,  in  order  to  perform  accurate  simulations,  aspects  of  both  premixed  and  non-premixed 
models  must  be  used.  With  this  motivation,  the  last  goal  of  part  I  was  to  validate  and  augment 
a  Combined  Conserved  Scalar /Level  Set  Flamelet  approach  based  on  a  multi-flamelet  model. 

In  the  following  sections,  the  methods  developed  to  meet  these  goals  will  be  described,  and 
the  results  of  the  simulations  performed  to  test  them  will  be  shown.  A  concluding  section  which 
summarizes  the  entire  project  will  follow  the  presentation  of  part  II  of  the  report. 

1.2  Non-Premixed  Turbulent  Combustion 

Reactions  in  non-premixed  combustion  in  technical  systems  most  often  are  limited  by  the  pro¬ 
cesses  that  mix  fuel  and  oxidizer.  While  macroscopic  structures  such  as  large-scale  vorticies  can 
force  some  of  this  mixing,  small-scale  mixing  is  required  to  achieve  maximum  heat  release  rates. 
This  requirement  exists  because  only  when  molecular  mixing  occurs  do  the  maximum  number 
of  reactant  molecules  come  into  contact  with  each  other.  The  mixture  fraction  variable,  Z, 
is  a  conserved  scalar  quantity  used  to  describe  this  mixing  between  fuel  and  oxidizer  streams. 
Thus,  an  appropriate  description  of  the  sub-filter  mixture  fraction  distribution  is  critical  for 
predictive  LES.  Normally,  the  statistical  moments  of  the  conserved  scalar,  namely  the  resolved 
scalar  and  the  sub-filter  variance,  are  used  to  obtain  a  PDF-based  description  of  the  sub-filter 
structure.  Treatment  of  the  resolved  scalar  with  a  transport  equation  is  straightforward,  but 
treatment  of  the  variance  is  problematic. 

Although  several  models  for  the  sub-filter  variance  exist  in  the  literature  [7,  8,  9],  none 
are  able  to  describe  the  sub-filter  structure  adequately.  Existing  models  rely  on  the  so-called 
local  equilibrium  assumption  that  equates  the  production  of  the  sub-filter  variance  to  the  dis¬ 
sipation.  In  turbulent  flows,  such  an  equilibrium  cannot  be  satisfied  instantaneously,  but  only 
in  an  average  sense.  In  spite  of  this  drawback,  the  dynamic  model,  which  makes  use  of  this 
assumption,  is  one  of  the  most  commonly  used  sub-filter  variance  models.  Adding  to  the  prob¬ 
lem,  the  numerical  implementation  of  this  model  involves  the  additional  ad-hoc  assumption  of 
directional  averaging,  which  is  required  to  remove  certain  singularities.  So  far,  no  convincing 
physical  explanation  has  been  provided  for  such  an  averaging  procedure.  A  final  difficulty  in 
correctly  describing  the  mixture  fraction  field  is  the  filter  size  itself.  Pope  pointed  out  that  LES 
is  an  incomplete  model  if  the  filter  size  can  be  arbitrarily  specified  [10].  In  this  case,  arbitrary 
specification  negates  the  possibility  of  setting  limits  on  the  requirements  for  a  sub-filter  variance 
model. 
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Work  performed  under  this  project  has  led  to  two  major  accomplishments  that  address 
these  issues.  The  first  is  the  development  of  a  new  method  for  describing  the  sub-filter  variance 
of  Z.  The  second  is  a  recursive  filter  refinement  (RFR)  procedure  designed  to  overcome  the 
arbitrary  nature  of  filter  specification. 

The  dynamic  modeling  (DM)  procedure  for  the  sub-filter  variance  mentioned  above  was 
proposed  by  Pierce  and  Moin  [8].  Based  on  the  assumption  of  local  equilibrium,  and  using  a 
modeled  dissipation  rate,  the  variance  model  can  be  written  as 

~Z?' ^  =  CA2(VZ)2 ,  (1) 

where  A  is  the  filter  width  and  C  is  a  coefficient  to  be  determined.  The  dynamic  modeling 
procedure  is  then  given  by 

fZ  -  ZZ  =  CA2 ( Vf  )2  -  CA^VZ)2  ,  (2) 

where  ~  denotes  test  filtering.  This  equation  is  solved  by  making  the  assumption  that  C  —  C, 
and  using  a  least  squares  estimation  technique. 

To  overcome  the  aforementioned  drawbacks  of  this  procedure,  a  new  method  has  been 
developed.  This  method  is  based  on  the  solution  of  a  transport  equation  for  the  second  moment 
of  the  mixture  fraction,  ( Z 2).  This  transport  equation  does  not  contain  any  production  terms, 
but  does  contain  the  scalar  dissipation  rate,  which  is  modeled  using 

X  =  ^(^-Z2),  (3) 

where  Cz  is  a  dynamically  determined  coefficient  and  r  is  a  suitable  turbulent  time-scale.  With 
the  solution  to  the  transport  equation 

^  +  V  -  (u Z2)  =  V  •  (ctVZ2)  +  V  ■  ( uZ 2  -  ^Z2)  -  ( 2a  (VZ)2  +  x)  ,  (4) 

the  sub-filter  variance  can  be  obtained  as 

Z*2  =  Z2-Z2.  (5) 

Unfortunately,  typical  dynamic  modeling  procedures  assume  azimuthal  homogeneity  and 
scale  independence  of  the  modeled  coefficient  ( Cz  in  Eq.  3).  To  overcome  this  issue,  the  dynamic 
localization  (DL)  procedure  is  used  to  formulate  an  integral  equation  that  then  is  solved  using  an 
iterative  technique  to  obtain  Cz  as  a  three-dimensional  field.  This  dynamic  modeling  procedure 
was  developed  assuming  local  equilibrium  and  tested  using  DNS,  as  well  as  LES  simulations. 
Results  from  shear-flow  calculations,  shown  in  Figure  1 ,  indicate  that  the  dynamic  model  over¬ 
predicts  the  scalar-to-mechanical  time-scale  ratio  compared  to  the  DL  model  by  a  factor  of  10. 
The  DL  model  consequently  predicts  a  higher  variance  across  the  shear-layer. 

The  second  accomplishment  in  this  part  of  the  program  was  the  development  of  the  recursive 
filter  refinement  (RFR)  method  [11].  In  this  method,  the  ratio  of  the  sub-filter  scalar  variance 
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Figure  1:  (Left)  Comparison  of  scalar-to-mechanica!  time-scale  ratios  R  computed  using  DM  (dashed)  and  DL 
(solid)  procedures.  (Right)  Comparison  of  variance  Z" 2  computed  using  DM  (dashed)  and  DM  (solid)  models. 


to  the  maximum  possible  variance  is  evaluated  and  used  to  describe  the  local  resolution  of  the 
mixture  fraction  field.  Since  the  maximum  local  variance  based  on  the  resolved  mixture  fraction 
is  Z(l-Z),  this  parameter  may  be  written 

Z”2 

a  =  Wz^)' 

When  a  is  high,  the  mixture  fraction  field  is  not  well  resolved,  and  a  large  demand  is  placed 
on  the  sub-grid  model.  When  a  is  low,  a  certain  level  of  field  resolution  is  guaranteed,  and  the 
importance  of  the  model  drops.  In  practice,  this  procedure  is  used  to  determine  the  distribution 
of  the  filter-size  that  is  suited  to  a  particular  flow.  The  distribution  of  a  can  be  evaluated  from 
an  instantaneous  LES  solution.  Everywhere  in  the  flow  domain  where  a  is  above  a  certain 
tolerance,  the  mesh  is  refined  in  order  to  guarantee  the  effectiveness  of  the  sub-grid  model. 

We  have  demonstrated  in  the  simulation  of  a  bluff-body  stabilized  flame,  which  was  exper¬ 
imentally  studied  by  Dally  et  al.  [12],  that  the  application  of  the  RFR  method  significantly 
improves  the  results  [11].  All  velocity  components,  scalar  and  temperature  profiles  predicted 
by  the  simulation  were  in  very  good  agreement  with  experimental  values.  Examples  of  temper¬ 
ature  and  carbon  monoxide  results  are  compared  with  experimental  data  in  Figs.  2  and  3. 

To  illustrate  the  complex  unsteady  fluid  dynamics  in  the  system,  a  ray-tracing  based  visu¬ 
alization  of  the  flame  is  shown  in  Figure  4.  It  is  clear  that  a  temporally  accurate  and  spatially 
well-resolved  simulation  is  essential  for  predicting  such  complex  flame  configurations. 


1.3  Premixed  Turbulent  Combustion 

Because  reactions  in  premixed  flows  occur  in  thin  fronts,  level  sets  have  been  suggested  as  means 
of  modeling  premixed  combustion  [13].  Level  sets  describe  iso-contours  of  field  variables.  They 
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Figure  2:  Radial  profiles  of  time-averaged  temperature  (solid  lines)  compared  with  experimental  data  (symbols) 
at  axial  locations  of  X  =  13,  30,  and  65  mm. 


Figure  3:  Radial  profiles  of  time-averaged  CO  mass  fraction  (solid  lines)  compared  with  experimental  data 
(symbols)  at  axial  locations  of  X  =  13,  30,  and  65  mm. 

are  governed  by  a  transport  equation  valid  at  the  iso-contour  of  interest  and  a  reinitialization 
procedure  that  sets  the  field  variable  away  from  the  iso-contour.  For  premixed  reacting  flows, 
the  so-called  G-equation  typically  is  used  to  describe  the  flame  front.  This  equation  differs  from 
a  simple  advection  equation  only  in  that  it  contains  a  source  term  describing  front  propagation. 

Since  the  inception  of  level  set-based  models,  a  number  of  issues  have  limited  their  appli¬ 
cation  to  LES.  The  first  such  issue  is  the  representation  of  the  burning  velocity.  To  describe  a 
variety  of  industrially  relevant  flames,  this  parameter  functionally  must  depend  on  local  quan¬ 
tities  such  as  the  equivalence  ratio,  pressure,  temperature,  and  turbulence.  A  second  issue  is 
filtering.  Because  the  level  set  transport  equation  is  valid  only  at  a  two-dimensional  surface  in 
the  three-dimensional  space,  volumetric  filtering  procedures  cannot  be  applied. 

In  meeting  the  goals  of  the  premixed  section  of  this  project,  a  new  filtering  procedure  for 
level  sets  was  developed  successfully.  That  procedure  led,  in  turn,  to  a  new  model  for  the 
turbulent  burning  velocity.  Finally,  based  on  simulations  performed  using  this  model,  a  new 
sensitivity  in  premixed  combustion  simulations  has  been  identified  as  an  area  in  which  addi¬ 
tional  work  could  improve  the  accuracy  of  these  types  of  simulations  greatly. 
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Figure  4:  Visualization  of  the  bluff-body  stabilized  flame  using  a  ray-tracing  technique. 


The  filtering  procedure  developed  for  the  G-equation  [14]  first  ’filters’  the  level  set  variable 
itself  by  integrating  it  over  a  computational  cell  volume  as  the  argument  of  a  Heaviside  function 

G(x,  t)  =  f  H(G(x  -  r,  t))dr.  (7) 

Jv 

The  Heaviside  function  separates  the  flow  in  two  regions  with  temperatures  below  and  above 
the  inner  layer  temperature.  The  filtered  field  then  is  an  average  over  these  regions.  Using  this 
new  variable,  a  filtered  transport  equation  may  be  derived  as 

r)G 

•  VG  =  sT,u  |VG|  .  (8) 

This  equation  has  two  unclosed  terms:  the  convection  velocity  vu  and  the  burning  velocity  st,u- 
The  convection  velocity  appearing  in  this  equation  is  an  average  conditioned  on  the  location  of 
the  flame  surface.  To  relate  vu  to  the  unconditionally  Favre-filtered  velocity,  v,  which  is  known 
from  the  solution  of  the  momentum  equations,  an  additional  model  was  developed.  Applying 
this  model  results  in  the  following  equation,  which  governs  the  filtered  flame  front  location 

+  v  ■  VG  =  =sT,u  |VG|  .  (9) 

While  standard  filtering  procedures  for  LES  do  not,  in  and  of  themselves,  introduce  errors  into 
the  equations  being  solved,  that  is  not  the  case  here.  The  need  to  model  the  convective  velocity 
used  in  the  filtered  level  set  equation  introduces  an  error  prior  to  numerical  solution.  It  can 
be  shown,  however,  that  this  error  is  proportional  to  some  power  of  the  filter  size  being  used, 


and  thus,  in  most  cases,  it  is  negligible.  Hence,  this  new  filtering  procedure  is  a  significant  result. 

Interestingly,  because  of  this  new  procedure,  a  propagation  term  proportional  to  the  cur¬ 
vature  of  the  mean  front,  which  appeared  in  earlier  formulations  of  the  filtered  G-equation,  is 
no  longer  present  in  Eq.  (9).  This  difference  is  important  because  the  curvature  term  has  a 
stabilizing  effect  on  the  flame  front  and  can,  therefore,  lead  to  a  reduction  in  resolved  flame 
wrinkling.  In  the  filtered  equation,  this  effect  instead  must  be  captured  by  the  model  for  the 
turbulent  burning  velocity,  the  second  unclosed  term. 

A  model  for  the  turbulent  burning  velocity  that  includes  this  stabilizing  effect  has  been 
developed  under  this  project  [14],  In  this  model,  the  front  propagation  velocity  is  described 
using  4  terms  as 

st,u  -  ( sl  +  st  -  Dk-  Dtk)  (10) 

which  account  for,  respectively,  the  resolved  laminar  propagation  speed  sl,  the  resolved  tur¬ 
bulence  induced  propagation  st ,  the  sub-filter  turbulent  contribution  Dk,  and  the  interaction 
of  resolved  curvature  and  sub-filter  mixing  Dtk.  Here,  D  is  the  diffusivity,  Dt  the  turbulent 
diffusivity,  and  k  the  resolved  curvature. 

For  completeness,  St  in  Eq.  (10)  is  developed  further  through  an  equation  describing  the 
magnitude  of  the  length  scale  separating  the  filtered  flame  front  location  from  the  exact  front 
location.  The  production  and  dissipation  terms  in  this  equation  are  assumed  to  be  in  balance, 
yielding  the  analytic  expression 


where  i/t  is  the  turbulent  viscosity,  Set  the  turbulent  Schmidt  number,  and  v'A  the  sub-grid  scale 
velocity. 

These  models  have  been  tested  in  simulations  of  the  premixed  turbulent  Bunsen  flames 
experimentally  studied  by  Chen  [15].  The  results  of  these  tests  are  shown  in  Figs.  5  and  6. 

Figure  6  shows  statistically  averaged  species  and  axial  velocity  results.  Species  such  as  H2O 
and  CO2  are  in  excellent  agreement  with  experimental  values.  Additionally,  only  very  small 
discrepancies  exist  in  the  axial  velocity  plots.  These  results  demonstrate  the  capabilities  of 
these  newly  developed  premixed  models. 

Figure  7  shows  data  from  similar  simulations  is  run  using  the  filtering  and  burning  velocity 
model  developed  here,  but  with  a  different,  supposedly  more  accurate,  heat  release  model.  In 
this  case,  the  transition  from  unburned  to  burned  gas  occurs  too  rapidly,  and  the  slopes  of  the 
species  plots  are  too  steep  near  the  flame  front.  These  results  indicate  that  there  is  more  work 
to  be  done  in  this  area.  Although  an  appropriate  filtering  procedure  and  propagation  model 
have  been  developed  for  level  set  based  methods  in  LES,  complete  confidence  in  predictions 
is  not  yet  possible.  Specifically,  the  simulations  run  under  this  project  have  demonstrated  the 
importance  of  the  models  describing  heat  release  around  the  flame  front,  which  need  to  be 
investigated  further. 
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Figure  5:  Contour  plots  of  F3  premixed  flame  simulation.  The  black  line  represents  the  flame  front.  With 
red  denoting  high  and  blue  denoting  low  values,  the  quantities  and  ranges  shown  are,  from  left  to  right:  (a) 
Axial  Velocity  [—7.5,37.5]  —  (b)  Mixture  Fraction  [0, 1]  (c)  Total  Enthalpy  [-1200,  -200]  (d)  Temperature 

[300,  2100]  K. 


Figure  6:  Radial  profiles  of  time-averaged  axial  velocity,  CO2  mass  fraction,  and  H2O  mass  fraction  from  F3 
premixed  flame  simulation.  Combustion  effects  around  the  front  are  captured  accurately. 


1.4  Partially-Premixed  Turbulent  Combustion 

In  many  technical  combustion  applications,  pure  non-premixed  or  premixed  combustion  modes 
may  not  be  valid,  and  partially  premixed  combustion  must  be  considered.  For  example,  in  the 
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Figure  7:  Radial  profiles  of  time-averaged  axial  velocity,  CO2  mass  fraction,  and  H2O  mass  fraction  from  F3 
premixed  flame  simulation.  Here,  combustion  effects  around  the  front  are  not  accurately  modeled. 

stabilization  region  of  a  lifted  non-premixed  flame,  fuel  and  oxidizer  have  time  to  mix  before 
reacting,  and  front-like  behavior  develops.  To  describe  such  situations,  a  Combined  Conserved 
Scalar /Level  Set  Flamelet  model,  in  which  both  premixed  and  non-premixed  combustion  modes 
are  considered,  has  been  proposed  [16,  17].  In  this  formulation,  a  level  set  method  is  used  to 
describe  partially  premixed  flame  propagation  and  to  distinguish  between  burned  and  unburned 
regions.  The  reacting  gases  in  the  post-flame  region  axe  then  described  by  a  laminar  diffusion 
flamelet  approach. 

In  the  initial  phase  of  this  work,  directed  toward  model  validation,  two  areas  in  which  this 
formulation  suffered  were  identified.  The  first  area  involved  how  flamelet  models,  which  typi¬ 
cally  had  been  considered  in  non-premixed  contexts,  should  be  expressed  near  flame  fronts.  The 
second  area  was  the  reliance  of  the  method  on  a  flamelet  formulation  for  density,  which  would 
not  be  valid  on  the  unburned  side  of  the  flame  front  where  preheating  occurs.  Since  similar 
issues  were  encountered  in  premixed  modeling,  a  level  set  formulation  was  used  as  a  starting 
point,  with  the  goal  of  formulating  solutions  valid  in  both  premixed  and  partially-premixed 
regimes. 

First,  a  new  method  for  computing  the  filtered  density  was  devised.  In  this  method,  the  den¬ 
sity  is  determined  directly  from  the  temperature,  for  which  a  transport  equation  is  solved.  This 
temperature  equation  is  to  be  solved  only  in  the  unburned  region.  The  burned  gas  temperature, 
determined  from  the  steady  flamelet  library,  is  enforced  through  a  source  term  as  a  boundary 
condition  at  the  flame  front.  The  flame  position  is  still  determined  from  the  G-equation  and 
then  imposed  on  the  temperature  equation.  It  would  be  incorrect  to  solve  only  a  temperature 
equation  in  numerical  computations  of  premixed  turbulent  combustion.  This  model  proved  to 
be  very  successful  and  was  used  in  the  premixed  computations  described  in  section  1.3. 

Second,  a  description  of  a  flamelet  model  valid  near  flame  fronts  was  developed.  In  this 
model,  a  PDF  describing  the  likelihood  of  finding  burned  gas  around  the  flame  front  is  presumed. 
This  PDF  is  Gaussian  in  shape,  with  the  mean  set  by  the  local  value  of  the  filtered  level  set 
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variable.  This  PDF  can  be  used  to  compute  the  temperature  equation  source  term,  but  more 
generally,  the  filtered  value  of  any  reactive  quantity  <j>  may  be  computed  as 

/0  roc 

P{G)dG  +  4>u  I  P(G)dG  —  4>bPb  +  <PuPu  >  (12) 

-oo  J0 

where  P(G)  is  the  PDF  and  subscripts  6  and  u  denote  burned  and  unburned  values,  respec¬ 
tively.  As  an  example  of  how  this  equation  is  applied,  consider  a  flamelet  model  that  presumes 
functional  dependencies  on  mixture  fraction,  Z,  and  total  enthalpy,  H .  Given  the  joint  PDF’s 
describing  those  dependencies,  a  filtered  quantity  conditioned  on  the  flame  front  can  be  written 
as 

4>  ( G ,  Z,H)=pb  JJ  <pb{Z,  H)P[Z ,  H)dZdH  +PuJJ  M%,  H)P{Z ,  H)dZdH  ,  (13) 

where  the  burned  and  unburned  values  of  the  scalar  are  computed  using  the  presumed  PDF 
describing  reactions  around  the  flame  front.  Because  of  the  success  of  this  formulation,  it  was 
used  with  slight  modifications  to  compute  species  concentrations  in  the  premixed  Bunsen  flame 
described  in  section  1.3.  Those  results  indicated  that  such  a  model  works  very  well  when  the 
PDF  describing  reactions  around  the  flame  front  is  correct. 

These  new  methods  were  applied  within  the  Combined  Conserved  Scalar/Level  Set  Flamelet 
Model,  which  was  tested  through  simulations  of  a  partially-premixed  case.  Specifically,  the  ex¬ 
periment  by  Spadaccini  et  al.  [18],  a  well  documented  gas-phase  dump  combustor,  was  used. 
This  experiment  mimics  a  simplified  aircraft  engine  combustor  where  fuel  and  air  enter  the 
combustion  chamber  separately,  but  the  flame  lifts  off  from  the  nozzle  and  is  stabilized  through 
multiple  recirculation  regions. 

Figure  8  shows  contour  plots  from  this  simulation.  These  figures  show  the  diffusive  nature 
of  the  burning  in  most  of  the  flame.  Without  the  indicated  level  set,  however,  predicting  the 
flame  lift-off  is  very  difficult.  Figure  9  shows  predicted  CO  mass  fractions  at  two  different 
cross-sections  in  the  combustor  compared  with  experimental  data.  Although  for  overall  lean 
combustion  CO  mass  fractions  are  very  sensitive  to  flow  conditions,  experimental  data  are 
matched  well. 

In  summary,  the  Combined  Conserved  Scalar /Level  Set  Flamelet  Model  has  been  demon¬ 
strated  to  be  successful  in  simulating  partially-premixed  combustion.  As  might  be  expected, 
however,  it  suffers  from  the  same  sensitivities  to  the  PDF  describing  the  flame  front  that  pre¬ 
mixed  formulations  do.  More  rigorously,  this  sensitivity  is  a  dependence  on  the  length  scale 
associated  with  the  PDF  in  Eq.  12.  With  further  development  of  the  methods  used  to  describe 
this  length  scale,  partially-premixed  simulations  will  become  only  more  accurate. 
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Figure  8:  Contour  plot  of  partially  premixed  dump  combustor  simulation.  The  black  line  represents  the  partially 
premixed  interpretation  of  the  flame  front.  With  red  denoting  high  and  blue  denoting  low  values,  the  quantities 
and  ranges  shown  are,  from  top  to  bottom:  (a)  Mixture  Fraction  [0, 1]  (b)  Temperature  [300, 2100]  K. 


Figure  9:  CO  mass  fractions  predicted  by  a  diffusion-only  and  a  combined  diffusion/level  set  model  at  two  axial 
locations  in  the  combustor. 

2  Development  of  Reduced  Kinetic  Mechanisms  for  JP-8  and 
Transportation  Surrogate  Fuels 

2.1  Introduction 

For  the  application  of  LES  in  combustion  simulations  of  real  technical  devices,  the  chemical 
details  of  the  fuel  have  to  be  considered.  Computations  with  detailed  chemistry  can  be  done, 
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for  instance,  for  natural  gas  but  become  very  complicated  for  the  class  of  fuels  typically  used  in 
transportation.  Gasoline,  diesel,  and  jet  fuels  are  refined  products  of  petroleum,  or  crude  oil, 
and  are  mixtures  of  hundreds  of  different  hydrocarbon  species  that  meet  a  certain  number  of 
general  physical  specifications  like  minimum  flash  point  or  maximum  freeze  point  temperatures, 
heat  of  combustion,  boiling  range,  and  volatility.  Thus,  there  can  be  significant  variations  in 
fuel  properties  and  compositions  from  various  feedstocks  and  refining  processes.  Representative 
samples  of  these  fuels,  called  reference  fuels,  are  used  in  the  industry  to  validate  engine  design 
and  performance.  However,  it  now  is  acknowledged  widely  that  much  simpler  surrogate  fuels 
should  be  used  in  laboratory  experiments  and  computational  simulations  to  improve  repro¬ 
ducibility  and  fundamental  understanding  of  the  combustion  processes. 

Surrogate  fuels  are  defined  as  mixtures  of  a  limited  number  of  hydrocarbon  species  that 
will  have  similar  physical  and/or  chemical  properties  to  those  of  the  real  fuels.  Among  the 
most  important  properties  that  might  have  to  be  matched  are  the  heat  release,  which  is  essen¬ 
tially  reflected  in  the  C/H  ratio,  the  density,  the  boiling  range,  the  tendency  to  form  pollutants 
such  as  NOx  and  soot,  burning  velocities,  and  auto-ignition  characteristics.  Most  of  the  hy¬ 
drocarbons  in  transportation  fuels  are  members  of  the  four  major  classes  of  hydrocarbons: 
paraffins,  aromatics,  naphthenes,  and  olefins.  Recently,  six  chemical  classes  have  been  isolated 
that  are  presumed  to  form  the  minimal  set  of  constituent  types  necessary  to  describe  the  chem¬ 
ical  composition  of  real  fuels  [19],  namely  iso-paraffins,  normal  paraffins,  single  ring  aromatics, 
cyc.lo-paraffins,  olefins,  and  multi-ring  aromatics. 

The  selection  of  components  for  a  surrogate  fuel  is  controlled  by  its  applications  and  ac¬ 
curacy  requirements.  The  approach  is  first  to  define  possible  candidates  for  surrogate  fuel 
ingredients.  These  candidates  should  be  as  close  to  the  real  fuel  components  as  possible.  How¬ 
ever,  the  actual  choice  is  restricted  to  those  fuels  for  which  either  detailed  chemical  schemes 
exist,  or  fuels  for  which  a  fair  amount  of  experimental  data  are  available,  and  detailed  mecha¬ 
nisms  can  be  constructed  with  some  confidence.  A  number  of  possible  candidates  belonging  to 
the  six  chemical  classes  presented  above  have  been  proposed  recently  [19].  Iso-octane  has  been 
proposed  as  an  iso-paraffin,  n-heptane,  n-hexadecane,  and  n-decane  as  normal  paraffins,  toluene 
and  different  xylenes  as  one-ring  aromatics,  methylcylcohexane  as  a  cyclo-paraffin,  1-pentene 
as  an  olefin,  and  a-methylnaphthalene  as  a  multi-ring  aromatic. 

Detailed  mechanisms  for  some  of  the  components  have  been  developed  by  Curran  et  al.  for 
iso-octane  [20]  and  n-heptane  [21],  by  Bikas  and  Peters  [22]  for  n-decane,  and  by  Pitsch  for  a- 
methylnaphthalene  [23].  Also,  semi-detailed  and  detailed  mechanisms  for  certain  mixtures  have 
been  provided.  For  example,  a  detailed  mechanism  for  mixtures  of  iso-octane  and  n-heptane, 
the  so-called  primary  reference  fuel  mixtures,  has  been  developed  by  Curran  et  al.  [24],  and 
semi-detailed  or  lumped  mechanisms  of  multi-component  jet  fuel  surrogates  have  been  proposed 
by  Violi  et  al.  [25]  and  Agosta  et  al.  [26]. 

Several  obstacles  have  to  be  overcome,  however,  to  use  these  mechanisms  as  fuel  surrogates 
in  computational  fluid  dynamics  (CFD)  simulations.  First,  in  constructing  a  multi-component 
surrogate  fuel,  several  detailed  mechanisms  of  individual  components  have  to  be  combined. 
Experience  shows  that  this  combination  is  a  formidable  task.  Another  complication  is  the  size 
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of  most  detailed  mechanisms.  Typically,  the  detailed  mechanisms  of  single  components  have 
hundreds,  or  even  nearly  one  thousand,  species.  A  detailed  mechanism  for  a  multi-component 
surrogate  will,  consequently,  also  have  perhaps  five  hundred  to  one  thousand  species.  For  most 
models  presently  used  in  LES,  the  maximum  possible  number  of  species  is  quite  restricted. 
Steady  flamelet.  models  can  tolerate  up  to  hundreds  of  species,  unsteady  flamelet  or  CMC  mod¬ 
els  perhaps  fifty,  and  transported  PDF  models  about  fifteen  to  twenty.  It  is  clear  that  chemical 
schemes  for  surrogate  fuels  cannot  be  used  unless  good  reduced  chemical  mechanisms  are  avail¬ 
able. 

The  top-down  approach  of  first  developing  detailed  mechanisms  for  individual  components, 
combining  these  mechanisms,  and  performing  the  reduction  is  very  time  consuming  and  inflex¬ 
ible.  For  example,  after  adding  a  new  component  to  the  detailed  mechanism,  the  reduction  has 
to  be  redone.  Also,  the  detailed  mechanism  has  to  include  all  reactions  for  auto-ignition  and 
formation  of  pollutants,  even  though  these  reactions  might  not  be  required  for  every  application. 

To  achieve  more  flexibility  in  the  surrogate  fuel  modeling,  we  have  proposed  and  demon¬ 
strated  a  different  approach.  Once  the  possible  components  of  different  surrogates  are  defined, 
validated  detailed  mechanisms  for  each  of  the  individual  components  are  identified.  These  de¬ 
tailed  mechanisms  are  reduced  independently  for  various  conditions  and  accuracy  requirements. 
The  skeletal  mechanisms  form  the  so-called  component  library.  Combining  the  corresponding 
elementary  skeletal  mechanisms  for  each  component  then  creates  a  variety  of  chemical  models 
for  surrogate  fuels.  The  skeletal  mechanism  for  the  multi-component  surrogate  then  is  validated 
through  a  comparison  with  experimental  data.  This  approach  makes  it  substantially  easier  to 
add  a  new  component  to  a  surrogate.  The  approach  also  greatly  simplifies  the  merging  of  the 
mechanisms  of  the  individual  components,  since  this  merging  now  is  done  at  the  level  of  the 
skeletal  mechanism. 

In  the  following  section,  the  method  first  will  be  reviewed  in  detail,  then  an  example  will  be 
provided,  showing  how  this  approach  can  be  used  to  design  a  skeletal  mechanism  for  a  surrogate 
fuel. 


2.2  Methodology 

As  mentioned  above,  detailed  kinetic  mechanisms  are  now  available  for  a  significant  number  of 
hydrocarbon  fuels.  However,  these  mechanisms  usually  are  designed  to  model  accurately  fuel 
oxidation  or  auto-ignition  over  a  large  domain  in  temperature,  pressure,  and  initial  concen¬ 
tration  space  and  can  include  thousands  of  reactions  and  up  to  one  thousand  species,  which 
prohibits  their  direct  implementation  in  flow  simulations.  Mechanism  reduction  aims  to  reduce 
those  mechanisms  to  the  smallest  possible  size  that  still  leads  to  an  accurate  prediction  of  a 
number  of  chosen  targets  over  the  domain  of  interest.  The  first  major  step  of  a  reduction 
process  identifies  all  species  and  reactions  that  have  little  or  no  influence  on  the  targets  and 
removes  them  definitively  from  the  mechanism.  The  number  of  differential  equations  that  have 
to  be  solved  is  reduced  by  reducing  the  number  of  species  and  the  computation  of  source  terms 
is  accelerated  by  removing  unimportant  reactions.  This  step  reduces  detailed  mechanisms  to 
the  so-called  skeletal  level. 
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As  most  of  the  computational  savings  come  from  the  reduction  of  the  number  of  differential 
equations  to  be  solved,  it  is  essential  to  be  able  to  identify  species  that  do  not  play  an  impor¬ 
tant  role  during  the  chemical  processes  correctly.  To  gain  insight  into  the  complex  chemical 
dynamics,  Bendtsen  et  al.  [27]  introduced  a  reaction  matrix  P  defined  at  any  time  t,  whose 
elements  P{j  correspond  to  the  production  of  species  j  from  all  reactions  involving  species  i  as 
a  reactant  at  time  t.  This  matrix  quantifies  the  interactions  existing  between  species  and  is 
used  to  select  important  species  iteratively.  The  selection  process  starts  with  a  set  of  major 
products  or  reactants.  Species  that  contribute  more  than  a  certain  percentage  to  the  overall 
removal  or  production  of  any  of  these  important  species  are  included  in  the  set.  This  procedure 
is  repeated  until  no  more  species  are  included.  The  results  were  used  to  generate  pathway  plots 
at  several  times  showing  graphically  the  conversion  of  fuel  into  main  intermediates  and  then 
into  main  products. 

Tham  et  al.  [28]  used  the  same  reaction  matrix  to  select  an  initial  pool  of  species:  for  each 
of  the  major  products  or  reactants,  an  additional  set  of  important  species  is  selected  by  going 
through  the  reaction  matrix  following  the  path  that  connects  one  species  to  the  next  that  is 
coupled  most  strongly  with  it.  A  second  step  selects,  for  each  of  those  species,  a  subset  of 
reactions  such  that  a  certain  percentage  of  the  total  formation  or  destruction  of  that  species  is 
kept  in  the  skeletal  mechanism.  Any  additional  species  appearing  in  these  reactions  then  are 
added  to  the  list  of  important  species,  and  this  last  step  is  repeated  until  no  more  species  are 
added.  Reactions  involving  large  heat  release  also  are  added  to  the  mechanism.  No  threshold 
is  needed  for  the  selection  of  the  initial  set  of  species.  This  method  is  efficient  in  selecting 
important  reactions  but  tends  to  retain  a  large  number  of  species  for  a  given  accuracy. 


Recently,  Lu  et  al.  [29]  used  a  species  selection  procedure  similar  to  the  one  proposed 
by  Bendtson  et  al.  [27]  but  applied  it  to  the  automatic  generation  of  skeletal  mechanisms 
with  various  complexities  and  domains  of  applicability.  Given  a  threshold  parameter  e,  a 
directed  relation  graph  (DRG)  is  constructed,  whose  nodes  correspond  to  species  present  in  the 
mechanism.  There  exists  a  directed  edge  from  species  A  to  species  B  only  if  the  normalized 
interaction  coefficient  tab  that  quantifies  the  dependence  of  species  A  on  species  B  is  greater 
than  t.  The  interaction  coefficient  is  defined  as 


tab  = 


£i=l,f  WA.iVjSBi} 
£»=i,/  WA,m\ 


(14) 


where  the  subscript  i  designates  the  2th  elementary  reaction,  2/4,1  the  stoichiometric  coefficient 
of  species  A  in  reaction  i,  u>i  the  production  rate  of  species  A  for  the  reaction  i,  and  5gl  is  the 
Kronecker  symbol 


.  _  J  1,  if  the  ith  reaction  involves  species  B 

1  \  0,  otherwise 

If  A  is  to  be  kept  in  the  skeletal  mechanism,  then  B  has  to  be  kept  too.  To  initiate  the  selec¬ 
tion  process,  the  user  must  specify  a  set  of  starting  species,  which  may  be  just  a  single  species, 
namely  the  fuel.  Species  kept  are  those  that  are  reachable  from  species  included  in  the  initial  set. 
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This  method  was  applied  to  auto-ignition  and  Perfectly  Stirred  Reactor  (PSR)  simulations. 
The  pre-defined  domain  of  applicability  of  the  skeletal  mechanism  to  be  generated  is  densely 
sampled  in  pressure,  equivalence  ratio,  temperature,  and  time.  The  DRG  selection  method 
is  applied  to  each  of  those  points,  and  the  final  resulting  skeletal  mechanism  is  the  union  of 
the  sets  of  species  obtained  at  each  sample  point.  The  selection  process  is  fast,  requires  a 
single  evaluation  of  the  solution  using  the  detailed  mechanism,  and  depends  on  one  unique 
user-defined  parameter.  However,  a  major  drawback  of  this  method  is  that  the  value  of  e  is  not 
directly  related  to  an  error  measure.  One  value  for  e  applied  to  different  mechanisms  will  return 
disparate  species  reduction  ratios  and  accuracies.  Moreover,  the  DRG  method  assumes  that 
every  species  selected  is  equally  important  and  that  the  set  of  strongly  coupled  species  associ¬ 
ated  with  each  important  species  has  to  be  kept  entirely,  which  may  not  be  necessary.  These 
issues  are  addressed  in  the  following  sections,  where  a  modified  version  of  the  DRG  method  is 
proposed  and  evaluated. 


2.2.1  Directed  Relation  Graph  Method  with  Error  Propagation  [30] 

For  each  species  A  present  in  a  kinetic  mechanism,  a  set  of  primary  dependent  species  can 
be  defined,  consisting  of  those  species  that  appear  explicitly  in  reactions  involving  A.  The 
strength  of  the  interaction  between  A  and  species  B  of  this  primary  dependent  set  is  defined  by 
the  interaction  coefficient  tab  defined  in  Eq.  14.  A  directed  relation  graph  can  be  constructed 
as  stated  by  Lu  et  al.  [29]  using  e  =  0.  If  species  B  is  not  in  the  primary  dependent  set  of  A, 
then  tab  —  0. 

Species  A  will  depend  not  only  on  the  species  contained  in  its  primary  dependent  set,  but 
also  on  all  species  directly  related  to  its  primary  dependent  set.  However,  species  C  interacting 
with  species  A  through  a  third  species  B  is  important  for  A  only  if  it  is  important  for  B  and 
if  B  is  important  for  A.  To  quantify  this  more  complex  coupling,  we  define  a  path-dependent 
coefficient  rAB,i  °u  path  i  from  A  to  B  as  being  the  product  of  each  primary  interaction  coef¬ 
ficient  encountered  on  the  path. 

In  Figure  10,  for  example,  if  path  #1  is  A  — > 

B  — >  D,  the  coupling  coefficient  between  A  and  D 
is 

D4D,i  =  tab  ■  rsD  ■ 

If  an  error  is  made  on  D,  it  propagates  through 
B  and,  the  larger  the  coupling  coefficient  between 
A  and  D  is,  the  larger  the  impact  this  error  will 
have  on  A.  Finally,  we  define  the  generalized  in¬ 
teraction  coefficient  of  species  A  with  species  B  as 
the  maximum  path-dependent  coefficient  between 
A  and  B 

Rab  =  max  {rAB,i}  . 

all  paths  l 

For  example,  in  Figure  10,  A  depends  on  C 


A 

i"ab/  \tac 

B  C 

I'AC ^  \l"BD 

C  D 


Figure  10:  Interaction  graph  between  four  species; 
coefficients  r  correspond  to  primary  interactions. 
Species  A  depends  on  species  C  and  D  through  its 
interaction  with  species  B 
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with  coefficient 


Rac  -  max  {rAB  ■  rBc,  rAc}  - 

Each  species  is  associated  with  a  subset  of  species  that  is  sorted  in  order  of  the  coupling 
strength  existing  between  them.  For  species  A,  a  subset  of  the  most  important  species  with 
respect  to  some  threshold  value  e  can  be  formed  by  keeping  in  the  set  only  the  species  i  for 
which 

Rab  >  e  - 

The  smaller  e  is,  the  more  species  are  included  in  the  set,  so  that  species  A  can  be  predicted 
more  accurately. 

The  skeletal  mechanism  corresponding  to  one  particular  value  of  e  then  is  generated  as 
follows.  First,  a  set  of  target  species  is  defined.  This  set  includes  the  species  whose  evolution 
should  be  described  accurately  by  the  skeletal  mechanism,  usually  the  major  reactants  and 
products,  NO  and  NO2  if  we  are  interested  in  NOx  formation,  and  any  additional  species  of 
special  interest.  By  defining  some  threshold  e,  a  subset  of  important  species  for  each  of  the 
targets  is  selected  automatically  following  the  procedure  described  above.  The  complete  set 
of  species  to  be  kept  is  the  union  of  all  these  subsets.  Any  elementary  reaction  containing  a 
species  that  has  not  been  selected  is  removed  from  the  mechanism. 

The  major  difference  between  this  new  selection  method  and  the  DRG  theory  is  that  it 
allows  a  finer  selection  of  species  since  it  focuses  on  how  the  error  made  on  any  single  species 
will  propagate  to  the  targets.  Selected  species  are  not  equally  important  anymore,  as  those 
species  that  are  farther  from  the  targets  are  comparatively  more  important  for  the  species  to 
which  they  are  related  than  those  directly  linked  to  the  targets.  As  a  consequence,  the  choice 
of  the  targets  must  be  made  carefully  since  only  the  targets  will  be  guaranteed  to  be  described 
accurately.  The  level  of  accuracy  of  any  other  species  will  be  determined  by  the  target  selection, 
and  the  agreement  between  the  skeletal  and  detailed  mechanism  may  be  quite  poor  for  those 
species.  This  discrepancy  is  not  a  problem  as  long  as  the  user  is  aware  of  the  conditions  for 
which  the  skeletal  mechanism  has  been  developed. 

An  illustration  of  this  species  selection  process  is  given  in  Figure  11.  Species  A  is  the 
target.  Species  B  and  D  contribute  5%  and  4%,  respectively,  to  the  production  of  target  A. 
Species  C  has  a  10%  effect  on  species  B.  Computing  the  global  coefficients  for  each  species  and 
rearranging,  even  if  C  had  the  largest  direct  interaction  coefficient,  it  only  has  a  small  effect  on 
the  prediction  of  A,  and,  therefore,  it  should  be  removed  first,  whereas  species  D  that  had  the 
smallest  coefficient  is  kept  in  the  reduced  scheme. 

The  DRGEP  method  identifies  unimportant  species  and  removes  them  from  the  skeletal 
mechanism.  The  method  insures  that  the  most  important  steps,  associated  with  the  most  im¬ 
portant  species  are  retained.  The  global  coefficients  do  not  distinguish  between  production  and 
consumption  of  the  species,  as  they  involve  the  absolute  values  of  the  production  and  consump¬ 
tion  rates.  However,  removing  all  the  consumption  reactions  for  a  species  can  sometimes  lead 
to  a  important  concentration  build-up  for  this  species,  introducing  an  error  more  important 
than  what  was  predicted  by  the  coefficient  itself.  An  additional  check  on  the  production  and 
consumption  of  each  species  is  made  during  the  selection  process  that  keeps  or  removes  extra 
species  so  as  to  avoid  such  build-up  cases. 
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A 


A 


Figure  11:  Example  of  species  selection  using  DRGEP 

Once  unimportant  species  are  identified,  not  every  reaction  step  involving  important  species 
is  important.  A  second  level  of  reduction  addresses  this  issue.  If  a  step  does  not  contribute 
significantly  to  the  production  or  consumption  rate  of  each  species  appearing  in  this  reaction, 
it  can  be  removed  safely,  provided  that  the  overall  production  and  consumption  rates  of  each  of 
those  species  remain  close  to  the  original  ones.  Assuming  that  the  first  step  keeps  all  necessary 
species  needed  to  reproduce  the  detailed  mechanism  results  to  a  given  accuracy,  removing  extra 
unimportant  reactions  should  leave  the  number  of  species  unchanged.  The  following  procedure 
thus  is  applied  to  each  species  selected  by  the  first  step.  Considering  forward  and  backward 
reactions  separately,  normalized  production  and  consumption  rates  of  a  species  i  with  respect 
to  a  reaction  j  are  defined  by: 

=  |  max(0,  Vjj)u;j\ 

13  Y,k  |max(0,r'ifcVfc|  ’ 

_  |  min(0,  Vifiwj \ 

13  £Jmin(0,i/ifcVfc|  ' 


For  each  species  i  that  has  been  retained,  both  the  normalized  production  and  consumption 
rates  for  each  irreversible  reaction  involving  that  species  are  sorted  by  increasing  value.  For  a 
given  threshold  6,  there  exists  an  index  J  such  that: 


g 

Eall  j  fij 


>  S  and 


gjU (Qjj) 

Eall  j  C>j 


>5. 


Only  reactions  from  1  to  J  are  included  in  the  subset  of  reactions.  The  final  subset  of  reactions 
is  the  union  of  the  subsets  obtained  for  each  species.  The  threshold  S  either  can  be  a  global 
threshold  or  can  depend  on  the  interactions  between  the  species  and  the  targets,  and  hence  on 
the  propagation  of  errors. 
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2.2.2  Automatic  Reduction  Procedure 

The  previous  section  explained  how  a  skeletal  mechanism  can  be  generated  for  a  given  c  and 
5.  This  section  addresses  issues  concerning  the  data  used  to  select  species  and  reactions,  how 
to  measure  the  error  between  skeletal  and  detailed  mechanisms,  and  the  iterative  procedure 
followed  to  obtain  the  values  of  e  and  S  that  give  the  best  reduction  possible.  This  section  also 
provides  an  overview  of  the  whole  reduction  procedure  from  the  perspective  of  a  user. 

Space  Sampling  Detailed  mechanisms  are  generally  valid  over  a  much  wider  domain  in  pres¬ 
sure,  temperature  and  concentration  than  needed  in  practical  applications.  A  skeletal  mech¬ 
anism  will  be  constructed  to  match  the  detailed  one  only  over  narrower  well-defined  regions. 
Several  points  that  are  characteristic  of  this  restricted  domain  are  selected  on  which  the  skeletal 
mechanism  will  be  required  to  satisfy  a  certain  accuracy.  If  those  points  are  well-chosen,  the 
reduced  mechanism  should  remain  valid  inside  the  region  mapped  by  them.  These  points  can 
include  auto-ignition  of  homogeneous  systems,  diffusion  or  premixed  flames  at  various  pressures, 
temperatures,  and  equivalence  ratios. 

Error  Evaluation  The  most  accurate  way  to  evaluate  the  error  introduced  in  the  skeletal 
mechanism  is  to  compute  the  solution  using  this  skeletal  mechanism  and  compare  it  to  the 
solution  obtained  using  the  detailed  mechanism.  At  each  sample  point,  depending  on  the  kind 
of  simulation  done,  an  error  is  computed  for  temperature,  ignition  delay  time,  laminar  burning 
velocity,  and  concentrations  of  the  targets,  scaled  to  be  between  0  and  1 .  The  global  errors  are 
the  maximum  errors  over  all  sample  points,  and  each  of  them  has  to  satisfy  the  corresponding 
user-defined  accuracy  requirement  for  the  skeletal  mechanism  to  be  valid. 

Iterative  Process  It  has  been  shown  that  this  method  removes  species  so  that  the  global 
error  increases  as  the  number  of  species  decreases  in  a  monotonic  and  relatively  smooth  way. 
Thus,  the  optimal  skeletal  mechanism  can  be  obtained  through  a  bisection  search.  Starting 
from  two  values  of  e,  the  first  one,  emax,  associated  with  an  unacceptable  error,  the  second  one, 
emim  with  an  error  within  the  tolerances,  a  new  e  is  chosen  equidistant  to  emin  and  emax  and 
the  corresponding  skeletal  mechanism  is  constructed  and  evaluated.  Depending  on  the  quality 
of  this  new  mechanism,  either  emjn  or  emax  is  replaced  by  this  new  value,  and  the  procedure  is 
repeated  until  the  best  acceptable  mechanism  is  found. 

Even  if  the  optimal  value  for  e  is  different  from  one  mechanism  to  another,  the  order  of 
magnitude  usually  remains  the  same.  The  initial  boundary  values  for  e  can  be  chosen  accord¬ 
ingly,  and  the  whole  procedure  converges  in  only  a  few  iterations.  Moreover,  evaluating  the 
solution  on  a  skeletal  mechanism  is  much  cheaper  than  on  the  detailed  one,  and  solutions  for 
large  values  of  e  are  seldom  evaluated  entirely,  as  the  skeletal  mechanism  is  rejected  as  soon  as 
the  solution  at  one  point  exceeds  the  tolerated  error. 

2.2.3  User  Input 

The  code  has  been  written  in  Perl  (Practical  Extraction  and  Report  Language)  and  is  interfaced 
with  the  chemistry  code  FlameMaster  [31].  Aside  the  kinetic  mechanism  and  the  thermody¬ 
namic  and  transport  data  for  the  species  appearing  in  the  mechanism,  the  user  must  provide 
an  input  file  that  includes  the  physical  configurations  for  which  the  skeletal  mechanism  must 
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be  valid,  the  target  parameters  for  each  case  and  the  error  tolerances.  An  example  of  input  file 
is  given  in  Figure  12.  A  skeletal  mechanism  that  satisfies  the  error  tolerances  for  each  target 
and  each  case  is  generated  automatically. 

2.3  Examples 

In  this  section  each  step  that  one  must  follow  to  generate  a  skeletal  model  for  a  surrogate 
fuel  will  be  presented.  The  detailed  kinetic  mechanisms  used  for  this  example  have  been  ob¬ 
tained  from  Lawrence  Livermore  National  Laboratories.  Taking  all  the  mechanisms  from  a 
single  source  presents  several  advantages.  First,  the  basic  C1-C4  mechanism  is  identical  for  all 
mechanisms,  which  greatly  simplifies  the  merging  of  the  skeletal  mechanisms  of  the  individual 
components.  This  merging  can  be  done,  in  general,  for  any  type  of  mechanism,  but  then  the 
validation  steps  become  much  more  crucial  and  should  be  done  with  extra  care.  Also,  the 
species  nomenclature  is  fairly  uniform,  which  avoids  any  duplicated  species  or  reactions  in  the 
final  combined  mechanism.  Finally,  the  thermodynamic  and  transport  data  are  similar  for  each 
detailed  mechanism.  Those  data  are  of  primary  importance  in  the  determination  of  reaction 
rates,  especially  for  reversible  reactions.  Using  one  set  of  thermodynamic  data  insures  that  the 
reaction  rates  in  the  combined  mechanism  will  remain  coherent.  Those  mechanisms  have  been 
validated  against  a  large  number  of  experimental  data.  Therefore,  we  can  have  confidence  in 
the  results  obtained  with  them. 

Although  taking  all  the  mechanisms  from  a  single  source  introduces  simplifications  in  the 
process,  there  is  a  major  drawback:  the  detailed  mechanisms  currently  available  do  not  include 
fuels  larger  than  iso-octane.  Jet  fuels  have  an  average  carbon  number  of  11,  and  appropriate 
surrogates  would  include  large  molecules  such  as  dodecane  or  butyl-cyclohexane.  On  the  other 
hand,  molecules  such  as  heptane,  iso-octane,  and  toluene  are  appropriate  for  gasoline  surrogates. 
A  surrogate  for  gasoline  was  designed  and  studied  experimentally  by  Gauthier  et  al.  [32].  Its 
composition  is  given  in  Table  1.  It  will  be  used  here  as  a  test  case  to  validate  the  method. 


vol% 

mol% 

Iso-octane 

63.0 

55.7 

Heptane 

17.0 

16.9 

Toluene 

20.0 

27.4 

Table  1:  Composition  of  Gasoline  Surrogate. 

The  detailed  mechanism  for  n-heptane  [21]  has  558  species,  and  the  detailed  mechanism  for 
iso-octane  [20]  has  850  species.  The  toluene  mechanism  has  been  extracted  from  a  mechanism 
for  HCCI  surrogate  fuel  by  Pitz  et  al.  [33]  and  consists  of  321  species. 

2.3.1  Validation  of  Detailed  Mechanisms 

Once  the  components  of  the  surrogate  are  defined  and  appropriate  detailed  mechanisms  have 
been  identified,  the  next  step  is  to  validate  those  detailed  mechanisms  against  experimental 
data.  Because  the  mechanisms  are  extremely  large,  with  hundreds  of  species,  the  only  validation 
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# - # 

#  GLOBAL  CHEMISTRY  PARAMETERS  # 

#  - # 

Common  { 

Mech:  . /gri.211.mech 

Thermo:  . /gr i.211.thermo.bin 

Fuels :  CH4 


Oxidizer:  02 
Global  Reaction: 

CH4  +  202  ==  C02  +  2H20 

> 


#  OD  CASES  # 

#  - # 

CommonOD  { 

Restol:  1.0e-9 

Abstol:  1.0e-9 

Noutputs :  100 
Pressure:  1.0e5 

> 

CaseOD  { 

Flame:  isochor 

Temperature:  950/100/1350 
Composition:  0.5  1  2 
Targets:  Tig/0.1  CH4/0/0.1 

> 

CaseOD  { 

Flame:  isobar 

Temperature:  950 
Composition:  1 

Targets:  CH4/0/0.1  C02/0/0.05 

> _ 

Figure  12: 


# - # 

#  ID  CASES  # 

#  - # 

CaselD  { 

Flame:  Counterflow  MixFrac 

TOxi:  400.0 

TFuel :  300.0 

VOxi:  -0.08745 

VFuel :  0.04821 

CompoOxi:  X->N2=0. 8;X->02=0.2; 

CompoFuel :  X->CH4=0 . 49 ; X->N2=0 . 51 

Pressure :  1 . 0e5 

Targets:  T/0/0.1  CH4/1/0.2  02/1/0.1 

> 


# - # 

#  REDUCTION  PARAMETERS  # 

#  - # 

NAME :  debug 

RECYCLE:  1  old_files 
COMPARE:  0 

SpecReduction  { 

Type:  standard 
Eval:  Linear  2  Stop 

> 

ReacReduction  { 

Type:  KeepAllSpecs 
Eval:  Bissection  0.9  1 

> 


Sample  input  file. 
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cases  that  can  be  done  at  this  point  are  zero-dimensional,  homogeneous  configurations  such  as 
shock  tube  (ST)  or  Plug-Flow  Reactor  (PFR)  experiments.  Figure  13(a)  shows  comparisons 
between  simulated  ignition  delay  times  and  experimental  data  from  Fieweger  et  al.  [34]  and 
from  Gauthier  et  al.  [32]  for  heptane.  The  simulation  results  match  the  experimental  data 
well,  but  the  mechanism  tends  to  over-predict  the  ignition  delay  times  at  high  temperatures. 
The  mechanism  for  toluene  had  to  be  modified  slightly  to  match  experimental  data  for  ignition 
times  and  species  concentrations.  Figure  13(b)  compares  the  simulated  time  evolution  of  major 
species  during  the  oxidation  of  toluene  in  a  PFR  with  experimental  data  from  Klotz  et  al.  [35]. 
The  agreement  here  is  very  good. 

2.3.2  Reduction  of  the  Detailed  Mechanisms  to  a  Skeletal  Level 

Each  detailed  mechanism  is  reduced  independently  to  a  skeletal  level  by  using  the  DRGEP 
method.  The  reduction  is  performed  for  shock  tube  and  PFR  simulations  over  a  wide  range 
of  initial  conditions,  and  the  maximum  error  tolerance  is  set  to  10%.  A  comparison  of  the 
number  of  species  before  and  after  the  reduction  is  shown  in  Table  2.  The  DRGEP  method  is 
very  efficient  at  removing  unimportant  parts  of  the  mechanism.  Figure  14  shows  a  comparison 
between  the  skeletal  mechanism  obtained  for  n-heptane  and  the  detailed  mechanism.  An  error 
of  10%  is  actually  reached  at  low  temperatures,  high  temperatures,  and  around  the  NTC  region. 


Component 

Detailed  Mechanism 

Skeletal  Mechanism 

Reduction  ratio 

Iso-octane 

850 

197 

4.3 

Heptane 

558 

170 

3.3 

Toluene 

321 

80 

4.0 

Table  2:  Reduction  of  Detailed  Mechanisms  for  the  Individual  Components. 


2.3.3  Combination  Step:  Skeletal  Mechanism  for  Gasoline  Surrogates 

A  skeletal  model  for  gasoline  can  be  obtained  by  combining  the  individual  skeletal  mechanisms. 
The  common  parts  of  each  mechanism  are  combined  and  the  parts  specific  to  the  individual 
components  are  added  automatically.  The  obtained  mechanism  for  gasoline  consists  of  352 
species.  Figure  15  shows  a  comparison  between  the  combined  mechanism  and  experimental 
data  by  Fieweger  et  al.  [34]  obtained  for  mixtures  of  heptane  and  iso-octane.  There  is  a  very 
good  agreement  between  the  two  mechanisms,  showing  that  cross-reactions  between  iso-octane 
and  heptane  have  little  influence  on  the  behavior  of  a  mixture.  Figure  16  compares  ignition 
delay  times  obtained  for  the  gasoline  surrogate  to  experimental  data  by  Gauthier  et  al.  [32]. 
Although  the  simulated  ignition  times  are  over-predicted  slightly,  the  overall  agreement  is  fairly 
good. 
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O  =  1.0 


(a)  Comparison  of  ignition  delay  times  between  LLNL  n-heptane  de¬ 
tailed  mechanism  (solid  line)  and  experimental  data  by  Fieweger  et 
al.  [34]  (open  symbols)  and  Gauthier  et  al.  [32]  (filled  symbols). 


Time  [ms] 


Time  [ms] 


(b)  Comparison  of  species  concentrations  in  PFR  between  LLNL  toluene 
detailed  mechanism  (solid  line)  and  experimental  data  by  Klotz  et 
al.  [35]  (symbols). 


Figure  13:  Validation  of  detailed  mechanisms  for  the  individual  components  of  a  surrogate. 
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1000/T  [1/K] 


(a)  Comparison  of  ignition  delay  times  of  skeletal  mechanism  (dashed  line) 
compared  with  detailed  mechanism  (solid  line)  and  . 


(b)  Absolute  value  of  error  of  skeletal  mechanism  compared  with  detailed 
mechanism. 


Figure  14:  Reduction  of  n-heptane  mechanism  -  Stoichiometric  ignition  delay  times  for  various  pressures:  p 
3.2  bar  (blue  line),  p  =  13.5  bar  (red  line)  and  p  =  42  bar  (green  line). 
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1000/T  [1/K] 


(b)  Pressure  =  40  bar 


Figure  15:  Reduction  of  primary  reference  fuels  mechanism  -  Comparison  of  ignition  delay  times  for  various 
octane  numbers:  ON  =  100  (red),  ON  =  90  (blue),  ON  =  80  (orange),  ON  =  60  (green),  ON  =  0  (magenta); 
detailed  mechanism  (black  solid  lines),  skeletal  mechanism  obtained  by  automatic  reduction  (colored  lined), 
experimental  data  (symbols).  94 


Figure  16:  Gasoline  ignition  delay  times  -  Comparison  between  simulations  (solid  line)  and  experimental  data 
for  gasoline  (filled  symbols)  and  gasoline  surrogate  (open  symbols)  at  different  pressures:  p  =  20  bar  (blue  line) 
and  p  =  55  bar  (red  line). 
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3  Conclusion 


In  the  first  part  of  this  work,  new  models  for  describing  sub-grid  quantities  in  reactive  LES 
settings  were  devised  and  validated.  These  models  included  a  new  model  for  the  sub-filter 
variance  of  a  conserved  scalar,  a  new  method  of  filtering  the  G-equation,  a  resolution  sensitive 
description  of  the  turbulent  burning  velocity,  and  a  flamelet  formulation  valid  near  premixed 
fronts.  The  models  have  been  shown  to  offer  improved  predictive  capability  through  application 
to  experimental  flames.  Furthermore,  the  new  models  have  contributed  to  the  identification  of 
a  significant  sensitivity  in  premixed  level  set  based  models.  Despite  the  development  of  these 
improved  models,  large  eddy  simulations  still  require  accurate  descriptions  of  chemistry  from 
outside  sources.  Most  often,  these  descriptions  come  from  the  solution  of  flamelet  equations, 
which,  in  turn,  require  accurate,  but  computationally  feasible,  chemical  mechanisms.  This  need 
motivates  the  second  part  of  the  work. 

In  that  second  part,  a  new  method  to  automatically  generate  skeletal  kinetic  mechanisms 
for  surrogate  fuels,  using  the  directed  relation  graph  method  with  error  propagation,  has  been 
proposed  and  evaluated.  These  mechanisms  are  guaranteed  to  match  results  obtained  using  de¬ 
tailed  chemistry  within  a  user-defined  accuracy  for  any  specified  target.  They  can  be  combined 
together  to  produce  adequate  chemical  models  for  surrogate  fuels.  A  library  containing  skele¬ 
tal  mechanisms  of  various  accuracies  and  domains  of  applicability  has  been  assembled.  New 
components  will  be  added  as  corresponding  validated  detailed  chemical  mechanisms  become 
available.  However,  the  resulting  combined  mechanisms  are  too  large  to  be  used  in  practical 
CFD  computation.  They  have  to  go  through  a  second  level  of  reduction  that  is  currently  under 
development.  The  strategies  used  in  this  second  stage  are  different  from  those  used  to  get  the 
skeletal  models  and  rely  on  automatic  lumping  of  species  and  the  use  of  Quasi-Steady  State 
Assumptions  to  reduce  the  number  of  differential  equations  that  have  to  be  solved. 

4  Participating  personnel 

During  the  reporting  period,  the  PI  and  two  graduate  students,  Perrine  Pepiot  and  Dirk  Veen- 
ema,  have  been  supported  within  the  present  research  grant. 

5  List  of  Submitted  and  Accepted  Publications  in  Reporting 
Period 

5.1  Refereed  Journal  Publications 

1.  Pitsch,  H.,  Large-Eddy  Simulation  of  Turbulent  Combustion,  Ann.  Rev.  Fluid  Mech., 
38,  pp.  453-483,  2006. 

2.  Pitsch,  H.,  A  Consistent  Level  Set  Formulation  for  Large-Eddy  Simulation  of  Premixed 
Turbulent  Combustion,  Comb.  Flame,  143,  pp.  587-598,  2005. 

3.  Raman,  V.,  Pitsch,  H.,  Large-Eddy  Simulation  of  a  Bluff-Body  Stabilized  Flame  Using  a 
Recursive-Refinement  Procedure,  Comb.  Flame,  142,  pp.  329-347,  2005. 
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6  Interactions  and  Transitions 

6.1  Selected  Presentations  at  Conferences  and  Seminars 

1.  Pitsch,  H.,  Large-Eddy  Simulation  of  Turbulent  Combustion,  Air  Force  Research  Labo¬ 
ratories,  October  2005 

2.  First  Workshop  on  Quality  Assessment  of  Unsteady  Methods  for  Turbulent  Combustion 
Prediction  and  Validation,  Frankfurt,  Germany,  “Validation  and  Verification  in  LES  of 
Reactive  Flows”,  June  2005 

3.  Pepiot,  P.,  Pitsch,  H.,  Systematic  Reduction  of  Large  Chemical  Mechanisms,  4th  Joint 
Meeting  of  the  U.S.  Sections  of  the  Combustion  Institute,  Philadelphia,  PA,  March  2005. 

4.  Workshop  on  Front  Propagation  and  Nonlinear  Stochastic  PDEs  for  Combustion  and 
other  Applications,  Montreal,  “Level  Set  Methods  in  Large-Eddy  Simulations  of  Premixed 
Turbulent  Combustion”,  January  2005 

6.2  Transitions 

In  the  last  year  we  have  worked  with  scientists  from  AFRL  to  make  different  computer  codes 
available  to  AFRL.  We  have  discussed  the  transition  of  the  software  used  for  automatic  re¬ 
duction  with  the  group  of  Dr.  Zelina.  Additionally,  we  are  in  the  process  of  transferring  the 
multi-physics  unstructured  LES  code  CDP  to  AFRL.  The  contact  people  at  AFRL  for  this 
effort  are  Dr.  Sekar  and  Dr.  Zelina. 

7  Patents  and  Inventions 

None 
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